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Rhythmic and sequential subdivision of the elongating vertebrate embryonic body axis into morphological 
somites is controlled by an oscillating multicellular genetic network termed the segmentation clock. This clock 
operates in the presomitic mesoderm (PSM), generating dynamic stripe patterns of oscillatory gene-expression 
across the field of PSM cells. How these spatial patterns, the clock's collective period, and the underlying 
cellular-level interactions are related is not understood. A theory encompassing temporal and spatial domains 
of local and collective aspects of the system is essential to tackle these questions. Our delayed coupling theory 
achieves this by representing the PSM as an array of phase oscillators, combining four key elements: a frequency 
profile of oscillators slowing across the PSM; coupling between neighboring oscillators; delay in coupling; and 
a moving boundary describing embryonic axis elongation. This theory predicts that the segmentation clock's 
collective period depends on delayed coupling. We derive an expression for pattern wavelength across the PSM 
and show how this can be used to fit dynamic wildtype gene-expression patterns, revealing the quantitative 
values of parameters controlling spatial and temporal organization of the oscillators in the system. Our theory 
can be used to analyze experimental perturbations, thereby identifying roles of genes involved in segmentation. 



INTRODUCTION 

During vertebrate development, segmentation of the con- 
tinually elongating embryonic body axis occurs rhythmically 
and sequentially from head to tail in a process termed somi- 
togenesis 11]. Somites are regularly sized cell clusters that 
bud off periodically from the anterior end of the posterior- 
most unsegmented tissue, the pre-somitic mesoderm (PSM), 
with a species-specific frequency. These transient, left-right 
symmetric structures are the embryonic precursors of adult 
bone and muscle segments, and defects in their formation 
lead to congenital birth defects Q. Underlying the morpho- 
genetic rhythm of somitogenesis, repeated waves of oscillat- 
ing gene expression sweep through the cells of the PSM from 
the posterior to the anterior O, see Fig.[T^ and Supplementary 
Movie 1 . These genetic oscillations are thought to slow down 
and arrest at different phases of their cycles at an anteriorly 
positioned arrest front that moves in concert with embryonic 
elongation |4 | (Fig. [TJ)), translating the temporal periodicity 
into a striped spatial pattern of gene expression. 

Given the existence of genetic oscillators in the cells of the 
PSM IS] 13, several questions still remain unanswered: how 
does a collective segmentation period arise from the popula- 
tion of individual oscillators, and what is the relation between 
the overall pattern of gene expression observed in the PSM 
and the oscillating expression at the cellular level. To investi- 
gate this, we develop a theoretical description based on phase 
oscillators and four key ingredients motivated by the biology: 
(i) a frequency profile along the PSM, slowing down the os- 
cillations, (ii) coupling of oscillators, (iii) a time delay in the 
information transfer between neighboring oscillators, and (iv) 
the existence of a moving front that arrests the oscillations at 



the anterior end of the PSM, while the posterior end moves 
due to embryonic outgrowth. This delayed coupling theory 
provides an excellent fit to the existing biological data, allows 
perturbations to the system to be analyzed in terms of underly- 
ing processes, and predicts how intercellular communication 
affects the collective period of the segmentation clock. Below, 
we introduce the elements of the delayed coupling theory. 



Phase oscillators 

To provide a simple picture of the segmentation process, 
Cooke and Zeeman proposed a clock and wavefront model 
more than thirty years ago fS]. However, to understand the 
role of collective processes in the emergence of dynamic gene 
expression patterns in the PSM a more detailed analysis is 
needed, for which methods from other pattern forming sys- 
tems can be boiTowed |9|. In particular, the periodic expres- 
sion of genes in the oscillating PSM cells can be described 
at tissue level using a set of phase oscillators, disregarding 
at this stage the underlying biochemical and genetic mecha- 
nisms that generate the oscillations and their pattern. In this 
phase description, each cell, or group of synchronous cells, 
is represented by an oscillator, and the state of each oscilla- 
tor is characterized only by its phase in the cycle of periodic 
gene expression. Oscillators with the same phase represent 
cells with equivalent expression level of cyclic genes. Pre- 
vious models described the PSM as a continuous oscillatory 
medium with a phase defined at each point of the PSM, see 
Supplementary data of |3| and |[T0l[IT]|T2][T3l. In this work 
we show that a phase description is sufficient to compute the 
overall spatiotemporal patterns of gene expression and the col- 
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FIG. 1: Representation of the pre-somitic mesoderm (PSM). Ante- 
rior is to the left and posterior to the right, (a) In situ hybridization 
1 5 1 showing the expression of deltaC mRNA in the zebrafish PSM 
(dorsal view), (b) Schematic PSM together with the already deter- 
mined segments — arrested — and the most recently formed pair of 
somites. The studied region lies between the arrest front and the pos- 
terior boundary, (c) Schematic representation of the signal gradient 
spanning the PSM (broken line). The frequency profile ui related to 
this gradient is depicted as solid purple line, using Eq. |2]l with the 
width (7 given in Table I. The length of the studied region is denoted 
by L. A linear array of A'^ coupled oscillators is indicated. 



profile, and hence its inclusion is purely phenomenological. 



Coupling of oscillators 

Recent theoretical works seeking to describe spatiotempo- 
ral patterns in somitogenesis using phase oscillators have not 
included coupling between oscillators llTOl [TTl [T2l [T3]| . al- 
though intercellular coupling has been considered in reduced 
models of regulatory circuits ll23l l24l 1251 . Here, coupling 
means that oscillators can influence the phase of their neigh- 
bors. Coupling is essential to stabilize tissue-scale patterns 
against the unavoidable noise present in biological systems 
Il25ll26ll27ll28l . and also to explain the re-synchronization of 
surgically inverted pieces of the PSM |4 |. Thus, in this work 
we propose a description based on coupled phase oscillators. 



Time delay 

Coupling between cells via signaling macromolecules, e.g. 
through the Notch pathway |l251 1251 W\ ISHl, involves syn- 
thesis and trafficking of such molecules within cells. These 
dynamics imply the existence of time delays, which have 
been recently estimated to be in the range of tens of min- 
utes in cell culture |f29l . Time delays in the coupling can 
have an impact on the self-organization of coupled oscillators 
Il23ll30l[3n[32l[33l . making their inclusion in our theory im- 
portant. For simplicity we will use here a deterministic time 
delay; a more realistic description would include a distribution 
of delays ll34l . 



lective period of the oscillations. 

Frequency profile 

It has been suggested that the arrest of the oscillations and 
the observed oscillating gene expression patterns are shaped 
by a spatial dependence of the frequency of the individual os- 
cillators 13 [ini [m [HI [H. A frequency profile could be 
controlled by the molecular gradients observed in the PSM, 
see Fig. [T];, e.g. the gradients of the growth factor FGF 
HI [H [151 [ISl or of Wnt signaling lUll [El. Several models 
have recently proposed regulatory mechanisms by which the 
genetic oscillations are affected by the gradients of signaling 
molecules along the PSM Ql] |20l [111 1211 ■ Motivated by the 
changing width of the stripes of gene expression in the PSM 
and the necessity that the oscillations slow down and finally 
stop at the arrest front, we include such a frequency profile in 
our theory. As with the assumption of cellular oscillators, our 
theory does not rely on the molecular origin of this frequency 



Moving borders due to embryonic elongation 

The embryo is a rapidly growing system, elongating about 
one somite length per oscillation cycle, which takes around 
25 minutes in zebrafish, 90 minutes in chick and 120 minutes 
in mouse. Cells are continuously in transit from the tailbud 
through the PSM, exiting it anteriorly as somites form. Addi- 
tionally, cell proliferation plays a role during elongation, but 
since in the PSM it has a stochastic character it can be con- 
sidered as a potential noise source l25l not otherwise signifi- 
cantly affecting the oscillatory dynamics |f35l, and we do not 
consider it further here. To correctly understand the forma- 
tion of patterns of gene expression and how the frequency is 
regulated, it is necessary to consider the geometry and bound- 
aries of the arena in which the process occurs. Here we ne- 
glect changes in the antero-posterior length of the PSM or 
the rate of axial growth, which occur during development at 
time scales larger than the somitogenesis period |[T3l[36l[37l . 
Consequently, both the arrest front and the posterior boundary 
move at the same velocity v, see Fig.[T]3. 
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(b) PSM REFERENCE FRAME 



FIG. 2: Two different coordinate systems to describe genetic oscil- 
lations in the PSM. (a) Lab reference frame: oscillators labeled by 
index i hold fixed positions Xi = ia, where a is the distance between 
oscillators. The PSM boundary moves posteriorly with velocity v as 
the embryo extends, (b) PSM reference frame: the PSM boundaries 
do not move but oscillators move through the PSM from posterior 
to anterior with velocity v. Oscillators are constantly relabeled us- 
ing symbols j to denote discrete positions relative to the arrest front. 
Dashed lines indicate the state of the system at slightly earlier times. 

RESULTS 

This section contains technical details of our theory. 
Readers who are more interested in the basic ideas and the 
biological justifications should note that a careful understand- 
ing of the equations is not a requisite to follow the arguments 
we introduce. 



Formulation of the delayed coupling theory 

Our model equations in the lab reference frame. Fig. |2^, 
consist of a lattice of discrete phase oscillators. This lat- 
tice comprises N oscillators in the antero-posterior direction, 
labeled by an index i. Each oscillator occupies a position 
Xi = ia along the PSM axis in the lab reference frame, where 
a is the characteristic distance between oscillators, i.e. the av- 
erage cell diameter. Hence the total physical length of the con- 
sidered system in the antero-posterior direction is L = Na. 

As the embryo elongates with velocity i;, the aiTest front is 
positioned at vt, where t is time. Oscillators anterior to the 
arrest front, Xi < vt, are arrested. New oscillators are added 
at the posterior boundary, situated at vt + L, as elongation 
proceeds. The oscillators in the studied region, with indices 
vt/a < i < {N + vt/a) are weakly coupled to their n nearest 
neighbors, denoted by the index k. In one dimension n = 2, 
while in a two-dimensional square lattice n — 4. The phase 
dynamics of the coupled oscillators can be described by 

e,{t) = + V sin Mt - n{t)) - e,it)] + at), 

k 

(1) 



where the dot denotes time derivative, di is the phase of os- 
cillator i, oji its intrinsic frequency, Si the coupling strength, 
Ti the time delay in the coupling and Q is a random variable 
with zero average representing different noise sources. Our 
objective in this work is to characterize the basic model in 
the synchronized state, and if not otherwise stated, we ignore 
the effects of noise. According to experimental evidence |27|, 
coupling strength is weak compared to other timescales in the 
system. This justifies the use of the sine function in the cou- 
pling, as it is the dominant term of any more general periodic 
coupling function |,38i . 

To specify the shape of the frequency profile we choose, for 
vt/a < i < {N + vt/a): 

u;,(i)=c^oo(l-e-(^'^-''*)/^), (2) 

where Uoo is a characteristic frequency scale of individual os- 
cillators, and (7 is a measure of the characteristic distance over 
which the frequency profile decreases from high to low values 
(see Fig.[T];). Below we will show that our choice of Eq. ^ is 
consistent with experimental observations, and we will deter- 
mine a from the width of the stripes of gene expression in the 
PSM. Qualitatively similar choices for the frequency profile 
have been used before ifTOl [TTI . For simplicity we have cho- 
sen the frequency to be strictly zero at the arrest front. Note 
however that this is not biologically necessary: a very low 
value of the frequency at the arrest front means a very large 
period of oscillation. If this period is much larger than any 
other time scale involved in the process, it determines in prac- 
tice an arrested oscillation, which can specify the downstream 
fixed pattern that eventually sets the position of somite bound- 
aries. Furthermore, since the oscillations can be coupled to a 
bistable system arising from opposing signaling gradients in 
the PSM, long period oscillations at the arrest front could be 
stopped by a bistable transition ETllSgl . 

For convenience we introduce the parameter lul = ^oo ( 1 — 
e~^/'^), that represents the intrinsic frequency of the oscil- 
lators at the posterior boundary of the system. Based on in 
situ experiments that show a largely uniform spatial expres- 
sion of cyclic genes in the tailbud at any given stage of the 
cycle (Fig. [T^), we define a uniform phase (Fig. [T]?) and ho- 
mogeneous frequency uj^ (Fig. [T];) in the tailbud, posterior 
to the notochord and the region we study here. This ho- 
mogeneity would be favored by the strong cell mixing |40|, 
and high, potentially saturating uniform levels of the signal- 
ing molecules that establish the gradient anterior to this region 
||4l[ll[T5l[l6l|T7l[T8|. Furthermore, the shape of lu, described 
by Eq. (j2| resembles the posterior branch of FGF receptor sat- 
uration proposed in |39|. 

The coupling could also be position dependent. In partic- 
ular, since the oscillators anterior to the arrest front stop cy- 
cling, they can not influence the active oscillators in the in- 
terval vt/a < i < {N + vt/a). We take this into account 
imposing = Sq = for i < vt/a. For simplicity, in this 
work we consider the coupling strength Si = e and the time 
delay = t to be constant posterior to the arrest front. 
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FIG. 3: Numerical simulation of segmentation using Eq. (jTJ in a 
growing two-dimensional geometry. Color intensity indicates the 
value of sin 6 of the phase 9: white is sin 9 = 1 and dark (red or blue) 
is sin 9 = —1. Vertical line indicates the position of the arrest front, 
with the oscillating PSM to its right (blue) and the arrested pattern to 
the left (red). Intrinsic frequency is a decaying function of the dis- 
tance to the black dot at the posterior boundary, causing the curvature 
of the stripes. We have used parameters determined for zebrafish at 
28°C, see Table I, and ujl = 0.224min^^, see main text for details. 
We have chosen to display three illustrative values of the time delay 
(a) T = min, (b) t = Tl/4 ^ 7 min and (c) r = 3Tl/4 ^ 21 
min. Delayed coupling affects the global frequency of oscillations 
according to Eq. The stripes of the oscillating PSM pattern and 
the segment length of the frozen pattern change accordingly. Movies 
available as supplementary material. 



Numerical simulations in two dimensions 

To gain insight into the role of delayed coupling in setting 
the period and the pattern of the genetic oscillations, as well 
as to illustrate the formation of realistic wave patterns within 
our theory, we performed computer simulations of Eq. ([T]i 
in a two dimensional geometry using different values of the 
time delay (Fig. [3] and Supplementary Movies 2, 3 and 4). 
Although the theory represents generic vertebrate segmenta- 
tion, here we use parameters from the zebrafish embryo. Ta- 
ble I. For the intrinsic frequency at the posterior we chose 
ujl = 0.224 min^^, corresponding to an intrinsic period of 
= l-K Iloj^ = 28 min. We show simulations with no time 
delay (r = min), a short delay compared to the intrinsic 
period (r = Ti/4 = 7 min), and a long delay close to the 
intrinsic period (r = ^TlI^ = 21 min). The latter is consis- 
tent with the experimental observation of tens of minutes for 
intercellular communication times ||29J . 

Fig. [3^ shows a snapshot of the simulation with no time de- 
lay. Not surprisingly, the collective period — the time needed 
to form one new arrested segment, and also the time after 



which the oscillating pattern in the PSM repeats itself — is 
unchanged with respect to the intrinsic period at the posterior 
The arrested segments have a length of 5 w 7 cell diameters. 
The case of short delay with respect to the period, t = Ti/4, 
qualitatively represents the situation in species with a rela- 
tively long segmentation peiiod, such as mouse. Fig.|3j) shows 
that in this case, the effect of the delay in coupling is to slow 
down the collective period (T = 39.1 min) with respect to the 
intrinsic period (T^ = 28 min). Further, the arrested segments 
are longer {S ~ 10 cell diameters) than in the case without de- 
lay, and there is a smaller number of gene expression stripes 
in the PSM, with increased size. Surprisingly, when the delay 
is made longer, the trends observed with the short delay are 
inverted. Fig.|3}; shows that for the time delay of r = 3Ti/4, 
the collective period (T = 23.5 min) is shorter than the intrin- 
sic period. Moreover, the arrested segments are also shorter 
(S* « 6 cell diameters) than in both previous cases, as are the 
stripes of oscillating gene expression in the PSM. 

This puzzling results show that delayed coupling introduces 
non-trivial effects to the system, large enough to be observable 
in real experiments. In order to understand these effects, in 
the following we perform an analysis of Eq. ([T]i, studying first 
the emergence of the collective peiiod from the parameters of 
the theory and then turning our attention to the spatial pattern. 
For this purpose, we will write Eq. ([T]i in a more convenient 
manner. 



PSM reference frame 

It is useful to consider the dynamics in the PSM reference 
frame, Fig.|2|3, where the oscillations can be characterized by 
a stationary phase profile and a collective frequency. For sim- 
plicity from here on we use a one dimensional description of 
the system, with n = 2. In the lab reference frame (Fig. |2^) 
the symbol i represents a fixed oscillator. In the PSM refer- 
ence frame (Fig. |2]3) we introduce the symbol j to label fixed 
discrete positions relative to the aiTest front. The label j runs 
from j = at the arrest front to j = N at the posterior bound- 
ary of the PSM. Discrete position j is occupied by different 
oscillators as the system evolves in time. For convenience we 
have included in the description the last arrested oscillator. 

In the PSM reference frame, the frequency profile is sta- 
tionary, LUj = ujooi^ ~ 6^^°/°^). Re-expressing Eq. (jlji in this 
PSM reference frame, an extra term describes the drift of the 
phase due to the movement of the cells relative to the PSM 
boundaries. The resulting phase dynamics are given by: 



(Pj{t) = UJj + V[lfij + l{t) - ipj{t)] + 

sm[ipk{t-T)- ipj{t)]. 



(3) 



2a2 



k=j+p±l 



Here ipj is the phase at position j relative to the arrest front 
and p ~ [vT /a] is the nearest integer tovr/a, representing the 
distance a cell moves during the time it takes for a signal from 
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FIG. 4: Global frequency Q, of somitogenesis as a function of time 
delay and coupling strength, (a) Dimensionless global frequency 
Q./u)L and time delay tujl are displayed for constant coupling e = 
0.07 (cell diameter)'^min~^, and intrinsic frequency ujl ~ 0.224 
min~^. Analytical solutions of Eq. ([sjl shown as blue lines: solid 
lines stable solutions of Eq. ([5j, dashed lines unstable solutions of 
Eq. l|5j. Blue dots correspond to numerical integration of the dis- 
crete model in two dimensions as given by Eq. l|TJ, for the two cases 
illustrated in Fig. [3j;,d. Red lines: global frequency as a function 
of delay in the continuum limit, showing its ranges of validity, (b) 
Global frequency as a function of time delay for different coupling 
strengths obtained from the solution of Eq. ^ with e = 0.11 (cell 
diameter) "^min^^ (green), e — 0.07 (cell diameter)^min^^ (blue) 
and £ = 0.03 (cell diameter) ^min~^ (red). Solid lines are stable so- 
lutions, dashed lines are unstable solutions. Dotted line at = 1 
corresponds to vanishing coupling, e = (cell diameter)^min^^. 



a neighbor to anive. Note that now the coupling is non-local: 
due to the time delay and cell movement, the neighbors of an 
oscillator with position j had positions j +p+l and j +p—l 
at the time the signal was sent. 



Steady state ansatz and collective frequency 

The oscillating gene expression pattern in the PSM repeats 
after a full period T = 27r/ri of oscillation ||3] |71, where 
is the collective frequency of the oscillation. This leads to the 
steady state ansatz (pj{t) — ilt + (j>j, where 0j is the station- 
ary phase profile describing the pattern in the PSM. With this 
ansatz we obtain from Eq. (j3]l: 

n = LUj+v[(t)j+i-(j)j] + ^ ^ sin [<j)k ~ <Pj - ^r] . (4) 

k—j+p±l 

The collective frequency of oscillations il is equivalent to the 
rate of somite formation. Note that the instantaneous frequen- 
cies di of individual oscillators depend on position and are in 
general different from il. 



Anterior boundary condition sets the segment length 

To determine the collective frequency we need to specify 
the boundary conditions, namely the conditions that (f>j ful- 
fills at the borders of the studied region (j = and j — N, 



Fig. [T]^). This boundary should not be confused with somite 
boundaries, which we do not discuss in this paper. 

At the arrest front {j = 0), the fact that ujq = and 
£o = implies with Eq. ^ that {(pi - (/)o) = fi/w = 2t:/vT. 
Thus, the anterior boundary condition determines the wave- 
length of the arrested pattern, which is the segment length 
S — 27r/((/)i — (po) = vT: the segment length is the distance 
advanced by the arrest front during one oscillation period |i8l. 



Coupling and delay affect the collective period 

At the posterior boundary of the PSM, we assume that new 
cells are added into the system with phase (f)]\i . To implement 
this we impose in Eq. Q p boundary conditions, (pj ~ 
for j = + 1, . . . , {N + p), accounting in this way for the 
effective non-locality of the coupling. We base this choice 
on the experimental observation of cyclic gene mRNA pat- 
terns, which maintain a smooth expression profile, and hence 
approximately homogeneous phase, across the interface be- 
tween tailbud and posterior PSM (e.g. Fig.[T^). 

Substituting the posterior boundary condition (t>N+i = 4>n 
in Eq. (|4]i we obtain a relation for the collective frequency of 
oscillations (see also |30, 31 , 32, 33|): 



= ojl — £ sin (rir) 



(5) 



The solutions to this equation are shown in Figs. |4^ and|4]5. 
Results from numerical simulations of Eq. ([T]) in two spatial 
dimensions show that the collective frequency indeed fulfills 
Eq. (|5]l, see blue dots in Fig.|4^. 

For a given set of parameters ujL,e, and r, Eq. (|5]l allows for 
multiple solutions for the collective frequency il. Independent 
measurement of coupling strength e, collective frequency lul 
and collective frequency il would allow the determination of 
possible values of the delay r consistent with Eq. (jSj. Exper- 
imentally, this can be done studying situations where the in- 
trinsic cellular oscillations are altered (modified ut^) or where 
the coupling strength is altered (modified e) and using the ob- 
served values of to fit r. 

A linear stability analysis following ll32l l33l reveals that 
when cos(rir) > the solution to Eq. (jsjl is stable, and un- 
stable otherwise, see continuous and dashed lines in Fig. |4] 
Consequently, multi-stability occurs for large values of e and 
r. As seen in Fig. |4^, for the biologically plausible parame- 
ters that we use in the figure, there is a small gap of time delay 
values between the first and the second stable branches of the 
solution. This happens around tlul — tt, which means that 
the delay is close to half the intrinsic period of the cellular os- 
cillators in the posterior PSM,t « Tl/2 = tt/ujl. For the pa- 
rameters in Fig.|4^, larger values of the delay always involve 
at least one stable solution. Note that values of the collective 
frequency equal to the intrinsic frequency, — uj^, aie only 
possible for delays equal to integer and semi-integer multiples 
of the intrinsic period T^: these solutions are stable in the case 
of integer multiples (r =integerxri) and unstable in the case 
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FIG. 5: Effects of noise in the delayed coupling theory. We include 
a white Gaussian noise as discussed in the text, (a) Delayed coupling 
is robust against the influence of noise. Parameters as in Fig.jSjj. (b) 
Impaired coupling results in segmentations defects. After initial syn- 
chronization with resulting segments not shown, coupling is turned 
off (e = (cell diameter)^min^^). The first segments have recog- 
nizable boundaries, but posterior segments are increasingly disrupted 
due to the effect of noise. Parameters as in (a). 



of semi-integer multiples (r = (integer+1/2) x TjJ). How- 
ever, stable solutions are possible for these latter delays, albeit 
with the collective frequency 17 different than the intrinsic fre- 
quency UJ^. 

Eq. (|5]l provides an explanation for the non-monotonic be- 
havior of the collective period observed in Fig. |3] Moreover, 
the simulation results coincide quantitatively with the predic- 
tion of Eq. (jSj), as shown by the three dots in Fig. |4^. Eq. (|5]l 
is biologically relevant: the collective frequency or period of 
somitogenesis emerges as a self-organized property, and de- 
pends not only on the intrinsic frequency of individual cells, 
but also on the coupling strength and the time delay (Fig. [4]i. 
Note that f2 does not depend on the specific shape of the fre- 
quency profile, and the period is set by the uniform phase cell 
population in the tail, which is the pacemaker of the whole 
oscillatory process. 
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FIG. 6: Phase profile in the PSM in the continuum limit, (a) Phase 
profile as a function of relative position, given by Eq. \\\\ . Left axis: 
phase relative to the arrest front. Right axis: corresponding num- 
ber of gene expression stripes. The green solid line corresponds to 
the set of parameters obtained from zebrafish data, see Table I, using 
Tt = 28 min and r = 21 min for illustration. Orange dotted line 
corresponds to ct = 6 cell diameters. (b,c) Waveform of the expres- 
sion pattern represented as sin i\>. (d) Number of stripes in the PSM as 
a function of a from Eq. (|7|, where a is the parameter describing the 
decay length of the frequency profile. Black solid line corresponds to 
parameters in Table I, with a variable. Green square dot: a obtained 
from zebrafish data, see Table I. Orange circular dot: mouse mode, 
see orange dotted curve in (a,c). Dotted and dashed curves corre- 
spond to higher and lower values of global frequency, which can po- 
tentially be affected by the intrinsic frequency, the coupling strength, 
or the time delay, see Eq. l[5](. (e) Zebrafish mode, reproducing panel 
(c) of Fig.|3]for comparison with panel (f). (f) Mouse/chick mode: 
zebrafish parameters as in (e), but with a sharper frequency profile, 
CT = 6 cell diameters. Only one wave of expression appears in the 
PSM, in contrast to the almost three waves in (e). Movies available 
as supplementary material. 



Delayed coupling keeps the oscillations synchronized 

We have seen how the presence of time delay in the cou- 
pling can have an important effect in the spatiotemporal pat- 
terns of gene expression. A critical biological function of 
intercellular coupling is to keep neighboring cells oscillat- 
ing in synchrony ll25l l26l l27l l28l . To demonstrate that de- 
lays in the coupling allow this function, and to and to show- 
case the role of noise in our theory, we simulate the phe- 
notype of a class of mutant embryos in which coupling is 
strongly reduced |[26l 1271 . To do this we include an addi- 
tive white Gaussian noise in the simulations, with zero mean 
and correlations {Ci{t)Ck{t')) = 2Q6{t — t')Sik, and choose 
^/2Q = 0.036 min~^ for illustrative purposes. With coupling 



as in Table I, the process is not disrupted by noise. Fig. [5^ 
and Supplementary Movie 5. When coupling is disrupted, 
the simulation exhibits posterior segmentation defects, Fig.|5]5 
and Supplementary Movie 6, resembling Delta-Notch mutant 
phenotypes in zebrafish [41] |42] |43] El HI). A more sub- 
tle feature of Fig. [S]? is the change in segment length after 
coupling has been disengaged. As soon as the coupling is re- 
moved, the effects of time delays are no longer present, and 
cells can oscillate at their intrinsic frequencies. Because at the 
time the coupling is turned off there is an established pattern 
in the oscillating PSM, it takes a few cycles for this informa- 
tion to be wiped out and to reach the new steady state value of 
the segment length. 

Although out of the scope of this work, our model provides 
a simple framework to study the effects of different kinds of 
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noise on segmentation. This interesting possibility remains 
open for future work. 



Spatial patterns of gene expression 

While the collective frequency describes the temporal 
regularity of somitogenesis, the spatial pattern of gene ex- 
pression in the PSM is characterized by the phase profile 
To evaluate the phase profile it is convenient to introduce a 
continuum limit where the spatial coordinate takes continuous 
values, denoted by x, replacing the discrete index j, see Meth- 
ods. The stationary phase profile see Fig. |6^, can be 
compared to quantitative experimental measurements of the 
pattern, such as the width of the stripes of gene expression 
reported in lfT2l . We define the wavelength A as the distance 
of two points in the PSM with a phase difference of 27r, see 
Fig.|7^. The wavelength is large close to the tail and becomes 
smaller close to the arrest front where it matches the segment 
length. Using the continuum formalism we find an expres- 
sion for the dependence of A with the position a; of the stripe's 
center relative to the arrest front 



CTlog 



sinh (A/2ct) 



7ri^-i(l + 77)-! + (A/2cr) e-^/"^ 



(6) 



Here 1/ and 77 are dimensionless parameters relating intrinsic 
frequency, coupling, time delay, elongation speed and the fre- 
quency profile, as defined in the Methods. In Fig.|7|3 we show 
the fit of Eq. (|6]l to the wavelengths obtained from the raw 
data in lfT2l : distances between consecutive points with equal 
level of herl expression in zebrafish embryos around the ten 
somite stage and raised at 28°C. The equation fits very well to 
the data, showing that our choice of Eq. Q for the frequency 
profile is consistent with observations. 



Parameter values 

From the fit to data obtained from wildtype zebrafish shown 
in Fig.|7j3 we determine L/a = 1.08 and + rj) — 57.8. 
We estimate the parameters L and a using the definitions of i/ 
and 7/ and the measured values of T and S, see Table I. Time 
delay affects both the collective frequency and the wavelength 
of the gene expression patterns. As we show in the Methods 
section, delayed coupling introduces a renormalization of both 
frequency and coupling strength. The effects of the time delay 
are thus included in the dimensionless renormalized parame- 
ters of Eq. (6), but the fit of spatial gene expression patterns 
does not allow the separation of the contribution of the time 
delay from that of the intrinsic frequency, and hence these two 
parameters remain undetermined from this fit. The intrinsic 
frequency at the posterior ujl, and the time delay t, are re- 
lated through Eq. (|5]l. Thus experimental determination of one 
would suffice to calculate the other if the coupling strength 
and collective frequency are known. 



T 


Period of somite formation 1371 


23.5 min 


n 


Global frequency, 2tt/T 


0.267 min"^ 


s 


Somite size (our own experimental estimation) 


6cd 


V 


Velocity of the arrest front, v = S/T 


0.255 cd/min 


e 


Coupling strength |27 | 


0.07 cd^min"^ 


L 


PSM length 


39 cd 


a 


Decay length of the frequency profile 


36 cd 


T 


Time delay 


Not determined 




Intrinsic freq. in the posterior PSM 



TABLE I: Parameters of the delayed coupling theory and their values 
in zebrafish embryo at 28°C near the ten somite stage. The first five 
parameters have been determined before or come from our obser- 
vations. PSM length L and decay length of the frequency profile a 
come from fits of our theory to experimental data in 1 12 |. The param- 
eters ujl and r could not be determined independently in this work, 
but they are related by Eq. l|5j. Note that parameters change with 
temperature and throughout development |37|. We choose as length 
unit one cell diameter (cd), in terms of which the distance between 
neighbor oscillators is a = led. 




FIG. 7: Wavelength of the pattern as a function of the position x in 
the PSM, where L is the length of the part of the PSM considered, (a) 
Schematic representation of the wavelength A. (b) Fit of Eq. ^ to 
the experimental data obtained from wildtype zebrafish in 1 1 2 1 . Best 
fit parameters are /i = 1.08 and + rj) — 57.8. (/i, u and rj are 
dimensionless quantities defined in the Methods.) 



From our estimated parameters in Table I the value of the 
frequency lul can be up to 30% higher or lower than the 
collective frequency fi, see Fig. |4^. For an intrinsic period 
Tl — 2tt/u)l around 28 min, this implies that changing cou- 
pling strength or delay time could situate the collective period 
in a range between 21 and 40 min, in qualitative agreement 
with the magnitude of period change from simulations of the 
genetic regulatory network model in |23 1 for two coupled cells 
B6l . Note that this period change is only possible due to the 
presence of delays in the coupling, both in our theory (see 
Eq. (|5])) and in the model in |46|. This difference in period 
is large and should be accessible to experimental observation, 
allowing at the same time for numerical determination of the 
values of the time delay r and the intrinsic frequency ujl. 



8 



DISCUSSION 



Methods we find: 



We have constructed a phenomenological theory describ- 
ing the tissue-level dynamics of the vertebrate segmentation 
clock employing phase oscillators to represent cyclic gene ex- 
pression in the cells of the PSM. As key ingredients of the the- 
ory, we considered (i) the existence of a frequency profile, (ii) 
coupling between oscillators, (iii) time delay in this coupling, 
and (iv) moving boundaries corresponding to embryonic elon- 
gation and the moving arrest front. Although these four ele- 
ments have been considered before, here we combine them 
in a unified framework. In this theory, tissue-level phenom- 
ena are generated by the interaction of cellular properties. For 
example, the collective frequency of oscillation of the PSM, 
related to the segmentation rate, depends on the intrinsic fre- 
quency at the posterior, the coupling strength and the time 
delay in the coupling, Eq. (j5]i; the spatial wavelength of gene 
expression stripes in addition depends on the shape of the fre- 
quency profile. Knowledge of the molecular underpinnings is 
not necessary for this mesoscopic description. By fitting the 
phase profiles obtained in our continuum limit to the existing 
data from a vertebrate embryo, we obtained a description of 
the tissue- and cellular-level processes controlling period and 
pattern in the system that is both quantitative and predictive. 
This framework can now be used to analyze experimental and 
evolutionary variants of embryonic segmentation or other per- 
mutations of growing, oscillating systems. 

Note that the basic relationship of a clock and wavefront 
type model for embryonic segmentation, as initially proposed 
by Cooke and Zeeman |8|, is that the length of a segment 
is the product of the arrest wavefront velocity and the period 
of the clock. In our description the population of oscillators 
create a pattern with a collective frequency, that together with 
the movement of the arrest front gives rise to a segment length 
consistent with the clock and wavefront picture. 



Variation of stripe patterns for different animal species 

We have compared our theory to zebrafish data, but it ap- 
plies equally well to other vertebrate species, since it does 
not involve species-specific details. The difference between 
what is termed a zebrafish mode of oscillation in somitogen- 
esis and a mouse/chick mode, observed also in medaka BtI . 
can be characterized as follows: in the zebrafish mode, sev- 
eral waves of gene expression sweep simultaneously through 
the PSM, i.e. multiple stripes of expression are detected in 
in situ experiments; in mouse/chick mode, only one wave is 
observed. The zebrafish mode applies also to snakes, where 
up to nine waves of gene expression have been observed lfT3l . 
Within our theory these different modes are characterized by 
the phase difference between the arrest front and the posterior 
border: the number of stripes of gene expression in the PSM 
is — 0(O))/27r, see Figs. |6^,b,c. From Eq. ( 11 1 in the 



Number of stripes 



vT {e^/<^-l)vT fis (eA'-l)s' 

(7) 

This expression can be written as a function of only two di- 
mensionless parameters: the ratio fi = L/a between the sys- 
tem length and the decay length of the frequency profile and 
the ratio s = S/L = vT/L between the segment length 
and the system length. The number of stripes is a decreas- 
ing function of both these ratios: smooth frequency profiles 
with long decay lengths, as well as small segment lengths, 
favor a large number of stripes of gene expression, as in the 
zebrafish mode, see Fig. |6}l. The coupling strength and time 
delay do not appear in Eq. (j7]i because we have neglected for 
simplicity higher order terms in e where they show up ex- 
plicitly. Note however that the collective period T = 2TT/il 
in Eq. (jTji does depend on both the coupling strength and the 
time delay through Eq. (jSj, see Fig. |6}i. Thus, in the same 
way that it may modify the collective period (as discussed in 
section ), the effect of delayed coupling can vary up to 30% 
the number of stripes of gene expression observed in the PSM 
compared to a system without coupling, see Fig. [3] 

A similar formalism for calculating the number of stripes 
has recently been published as Supplementary Material in 
IIT3I . The underlying theory was previously proposed in 
ifTOlfTTl . and is the same as our continuum theory, but without 
coupling, and hence without the effects caused by the delay in 
the coupling. However, in [| 13 1 no explicit choice for the shape 
of the frequency profile is made, hence the resulting formula 
for the number of stripes is a function of an unknown integral, 
rather than a closed formula. Our Eq. (|7]l allows for direct 
quantitative comparison with data. The choice of Eq. Q for 
the frequency profile comes from phenomenological observa- 
tions and it is not derived from the underlying molecular inter- 
actions of the signalling gradients with PSM cells. Neverthe- 
less, our function for the frequency profile is well supported 
by experimental data, see Fig. |7] 

It is important to note that a switch between modes can 
be achieved while preserving the timing of somitogenesis by 
changing the shape of the frequency profile: in Fig. ^ we 
show results of simulations of a mouse mode in zebrafish with 
all parameters given as in Table I except for cr, which is 6 
cell diameters instead of 36 cell diameters, see also Supple- 
mentary Movie 7. This implies that the number of stripes 
can change by changing the shape of the frequency profile 
while leaving the collective period and segment length un- 
affected. Previous hypotheses for the different modes in- 
clude changes in period, loss of stripe specific cyclic gene en- 
hancers, changes to the stabiUty of cyclic mRNA or different 
elongation velocities fT3l l47l l48l |49|. The delayed coupling 
theory indicates that changes to the frequency profile, poten- 
tially through changes to FGF or Wnt signaling gradients in 
the PSM, and different sizes of the PSM must be considered 
as well. This is consistent with recent experiments reported in 
IfTSl . where extra stripes of gene expression appear in a mutant 
with an expanded PSM. 
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Relation to regulatory network models 

Current regulatory network models for the genetic oscilla- 
tions in somitogenesis |[l9]|23]|50l|5l]|52]|53l, undergo a Hopf 
bifurcation — a generic mechanism by which oscillations can 
appear in a dynamical system — when varying some param- 
eters of the models, as for instance the transcriptional delays 
Il54l|55]|56l|57]|58l. Although it is also vaUd in more general 
settings, our Eq. ([T]l can be obtained as the phase equation as- 
sociated to the normal form of a Hopf bifurcation when varia- 
tions in the amplitude of the oscillations can be neglected, and 
as such it can in principle be derived from any of the dynami- 
cal systems associated with these regulatory networks follow- 
ing standard procedures 1381 l59l |60 1 . Hence our formulation 
represents a simplification that captures general features and 
properties of more detailed models. 

The mechanism arresting the oscillations at the arrest front 
is a different problem not addressed in our present work. 
While the above mentioned models undergo a Hopf bifurca- 
tion when varying one of their parameters, something com- 
pletely different (another kind of bifurcation triggered by the 
variation of a different parameter of the models, for instance) 
may be happening at the arrest front. The possibility that the 
oscillations are coupled to a bistable switch related to the sig- 
naling gradients in the PSM has been proposed 1211 |39l . In 
this scenario the arrest of the oscillators would not be a result 
of the intrinsic mechanism of the oscillations, but would result 
from an external signal. 

Implications of multistability 

Only stable solutions of our theory can be biologically rel- 
evant. In addition, we hypothesize that unique solutions are 
required to guarantee a robust behavior in the developing em- 
bryo. In the presence of multiple stable solutions for the 
collective frequency, fluctuations could drive the system to 
switch between these different states, with dramatic conse- 
quences for healthy development. For this reason, we con- 
jecture that if several time delays are consistent with a fit to 
experimental data, those yielding a unique value of the col- 
lective frequency should be favored. Biochemical evidence 
indicates that coupling time delays should be relatively short 
compared to other signaling processes in the vertebrate seg- 
mentation clock 1 29 1, thus likely precluding the observation 
of multistability in such an embryonic system. Multistabil- 
ity has been observed in other systems where coupling delays 
can be large, with applications in biochemistry [61], chem- 
istry |62]|63l, control theory ll64l or laser physics ||65[|661 . for 
example. 

Applications in somitogenesis and comparison to experiments 

Key quantitative experiments in vertebrate segmentation in- 
clude determination of segmentation rates [37J. and the anal- 



ysis of expression patterns from in situ experiments lfT2ll and 
fluorescent reporter genes |7 1. Our theoretical description al- 
lows for quantitative analysis of these experiments. 

The comparison to experimentally observed dynamic pat- 
terns of gene expression permits the determination of the 
model parameters, which are provided for wildtype zebrafish 
in Table I. Future studies in mutant embryos or embryos 
treated with different inhibitors will reveal which parame- 
ters are affected. The parameters in our model can be re- 
lated to different cellular functions such as molecular synthe- 
sis and trafficking of signals in the cell (coupling delay r); 
the strength of intercellular signaling (coupling strength e); 
the speed of a cell autonomous oscillator (intrinsic frequency 
ujl); changes in the signaling gradients responsible for the fre- 
quency profile (decay length a); and changes in the position 
of the arrest front (reflected by the system length L). Thus, 
analysis of experimental results using our theory can provide 
a deeper understanding of how molecular changes lead to new 
phenotypes from the altered collective dynamics of tissues. 

Our framework can be extended to other developmental 
processes that combine growth with a molecular clock. These 
are for instance fore-limb autopod outgrowth and patterning 
ll67l . or segmentation in short germ band insects, spiders, cen- 
tipedes, and other invertebrates that might form segments by 
a mechanism similar to the one we described ll68i i69l . 

Summary 

The delayed coupling theory describes spatiotemporal 
patterns of gene expression during morphogenesis in agree- 
ment with experimental observations. Most importantly, 
our phenomenological theory provides a unified quantitative 
framework relating the segmentation period and cyclic 
patterns of gene expression to underlying properties, such 
as the characteristics of intercellular communication, cell 
autonomous oscillations, the spatial profile of the slowing of 
the oscillators in the PSM, the rate of axial growth and the 
size of the PSM. Our results indicate that the specific spatial 
pattern of cyclic gene expression in the PSM does not affect 
the overall timing of somitogenesis, but intercellular commu- 
nication should be considered as a fundamental mechanism 
in setting the collective frequency of the segmentation clock. 



METHODS 
Continuum limit 

Starting from Eq. Q a continuum limit describing the evo- 
lution of the phase can systematically be derived for any value 
of the time delay. This continuum limit is valid when the typ- 
ical length scale of the modulations of the pattern is much 
larger than the distance between oscillators, a. The limit is 
obtained by letting the distance a tend to zero, while the total 
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number of oscillators N tends to infinite, in such a way that 
the length of the PSM, L = Na, remains finite and constant. 
In the continum limit, we require a finite coupling strength 

= lima^o s/a^ to exist, which implies that e scales as a^. 

The description based on discrete oscillators with phase 
(fij (t) at a distance aj from the arrest front (where j is a dis- 
crete label) is substituted by a description defined in a con- 
tinuous field spanning from x = to a; = L, where a; is a 
real positive value giving the distance to the arrest front of a 
point of the field with phase (p{x, t). The resulting continuum 
equation reads: 



ip{x, t) = uj{x) + v\Jip{x, t) + -^\l'^ip{x, t), 



(8) 



where v is the velocity of the arrest front, V denotes spa- 
tial derivatives (V = [d/dx) in one dimension, V = 
{d/dx, d/dy) in two dimensions, and son on), u;{x) is a posi- 
tion dependent effective frequency given by 



U>{x) — Lu{x) 



1 + 2Tnnec/ujL 



1 



(9) 



and ec = ec(l + 27rmec/wL)/(l + ^ct) is the effective cou- 
pling strength. The effect of the time delay appears through 
T and m = [tujl/2tt], the nearest integer to tllJl/2tt. In 
analogy with uij in the discrete case, the intrinsic frequency 



is defined as lu(x) 



.(1 



-X I a 



). Note that for simplic- 



ity we have assumed that the intrinsic coupling ec is constant 
throughout the PSM (as we did with e); it is straightforward 
to include a positional dependence by substituting by ec(a;) 
in all the previous expressions. 

We can simplify Eq. ([8]l using the steady state ansatz 
^{x, = fit + (i>{x) as we did in the discrete case: 



^ = uj{x) + vSI(\){x) + ^V^0(a;). 



(10) 

x)\,j,=Q = 



The boundary conditions for Eq. (lOi are 
and '^'l>{x)\x=L — 0. As in the discrete case, we assume 
that the phase is defined and uniform in the tailbud, 0(x > 
L) — (j>{L). This implies that at x — L all the derivatives in 
Eq. (fTo| vanish and fl = u}{L). In fact the right hand side 
of Eq. (|9|l coincides with the expression for fl obtained from 
solving Eq. Q after linearization around values of the delay 
r = i-Kin/uji^. In Fig|4^ we show in red the dependence of 
— u}{L) with T given by Eq. (j9]l for several values of to; 
note that almost the whole range of stable solutions of Eq. (|5]l 
(solid blue) can be well approximated by the continuum limit 
(red). 

Eq. ^ with u}{x) given by Eq. (j9| is valid when t w 
21X1X1 1 u>L for an integer m. A different equation can be ob- 
tained in the cases where t w 2'K(m + \/2)/ujl: it corre- 
sponds to the continuum approximation of the unstable solu- 
tions of Eq. (|5]l shown by the broken blue lines in Fig|4^. 

Eq. ( [TO) l can be solved and the corresponding phase profile 
reads: 



0(O = ^(l-r;)-i{(l-r;2) 



(11) 



where we have defined the dimensionless coordinate = x/L 
and parameters fi — L/a, v = Hi^a jv, and 77 = ec/2av; 
(/)(0) has been set to (/)(0) = to fix an arbitrary constant. 
Fig. [6^ shows the shape of this phase profile. 

The wavelength of the patterns of gene expression can be 
measured as a function of the relative position ^. In [12] this 
is done experimentally, using a definition of the wavelength A 
that in our notation can be expressed as the condition + 
A/2) — — A/2) = 2tt. In the limit of small coupling 
strength rj <C 1, we obtain a simple relation between the local 
wavelength of the pattern A and the position ^ along the PSM 
given by Eq. (j6]l of the main text. 
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